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, In this brief review we summarize a number of recent developments in the study of 

% I ■ vortices in Bose-Einstein condensates, a topic of considerable theoretical and experi- 

' mental interest in the past few years. We examine the generation of vortices by means of 

phase imprinting, as well as via dynamical instabilities. Their stability is subsequently 
examined in the presence of purely magnetic trapping, and in the combined presence of 
magnetic and optical trapping. We then study pairs of vortices and their interactions, il- 
lustrating a reduced description in terms of ordinary differential equations for the vortex 
centers. In the realm of two vortices we also consider the existence of stable dipole clus- 
ters for two-component condensates. Last but not least, we discuss mesoscopic patterns 
formed by vortices, the so-called vortex lattices and analyze some of their intriguing 
dynamical features. A number of interesting future directions are highlighted. 
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1. Abbreviations 

o BEC: Bose-Einstein condensate o 

o GP: Gross-Pitaevskii (Equation) o 

o MD: Molecular dynamics o 

o MT: Magnetic trap o 



OL: Optical lattice 
PR: Parrincllo-Rahman 
TF: Thomas-Fermi 

NLS: Nonlinear Schrodinger (Equation) 



2. Overview 

Bose-Einstein condensation (BEC) was theoretically predicted by Bose and Einstein 
in 1924 1,2,3a . It consists of the macroscopic occupation of the ground state of a 
gas of bosons, below a critical transition temperature T c , i.e., a quantum phase 
transition. However, this prediction was only experimentally verified after 70 years, 
by an amazing series of experiments in 1995 in dilute atomic vapors, namely of 
rubidium 4 and sodium 5 . In the same year, first signatures of the occurrence of 
BEC were also reported in vapors of lithium 6 (and were later more systematically 
confirmed). The ability to controllably cool alkali atoms (currently over 35 groups 
around the world can routinely produce BECs) at sufficiently low temperatures and 
confine them via a combination of magnetic and optical techniques (for a review 
see e.g., Ref. [7]), has been instrumental in this major feat whose significance has 
already been acknowledged through the 2001 Nobel prize in Physics. 

This development is of particular interest also from a theoretical/mathematical 
standpoint. On the one hand, there is a detailed experimental control over the 
produced BECs. On the other, equally importantly, there is a very good model 
partial differential equation (PDE) that can describe, at the mean field level, the 
behavior of the condensates. This model (which, at heart, approximates a quantum 
many-body interaction with a classical, but nonlinear self-interaction) is the well- 
known Gross-Pitaevskii equation 8,9 , a variant of the famous Nonlinear Schrodinger 
equation (NLS) 10 that reads: 

^ t = -|^A* + 5 |*| 2 * + K x t(r)*, (1) 

where * is the mean-field condensate wavefunction (the atom density is n — 
\^f(x, t)\ 2 ), A is the Laplacian, m is the atomic mass, and the nonlinearity coef- 
ficient g (arising from the interatomic interactions) is proportional to the atomic 
scattering length 7 . This coefficient is positive (e.g., for rubidium and sodium) or 
negative (e.g., for lithium) for repulsive or attractive interatomic interactions re- 
spectively, corresponding to defocusing or focusing cubic (Kerr) nonlinearities in 
the context of nonlinear optics 10 . Notice, however, that experimental "wizardry" 
can even manipulate the scattering length (and, thus, the nonlinearity coefficient g) 
using the so-called Feshbach resonances 11 to achieve any positive or negative value 



a We should mention that part of this overview section has significant overlap with the earlier 
review of two of the present authors on "Pattern Forming Dynamical Instabilities of Bose-Einstein 
Condensates" (Ref. [51] herein); it is included, however, here for reasons of completeness. 
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of the scattering length at will (i.e., nonlinearity strength in Eq. (1)). Moreover, 
the external potential V cxt can assume different forms. For the "standard" mag- 
netic trap (MT) usually implemented to confine the condensate, this potential has 
a typical harmonic form: 

Vmt = \m (lo 2 x x 2 + cu 2 y y 2 + cu 2 z z 2 ) , (2) 

where, in general, the trap frequencies oj x ,oj y ,uj z along the three directions are dif- 
ferent. As a result, in recent experiments, the shape of the trap (and, hence, the 
form of the condensate itself) can range from isotropic to the so-called cigar-shaped 
traps (see, e.g., Ref. [7]), to quasi- two-dimensional 12 , or even quasi-one dimensional 
forms 13 . Moreover, linear ramps of (gravitational) potential V ext = mgz have also 
been experimentally used 14 . Another prominent example of an experimentally fea- 
sible potential is imposed by a pair of laser beams forming a standing wave which 
generates a periodic optical potential, the so-called optical lattice (OL) 15 ' 16 ' 17 ' 18 , 
of the form: 



Vol = V 



,2 



/2irx , \ 2 {2ny \ 2 / 2nz 
L /A I 4- r-ns 2 -4- rh.. 4- rns^ 



COS y— r (px J ~r COS y— 1 <p y j + COS 



(3) 



where X x .y,z — -^lasor sm(6 X! y tZ /2)/2, Ai asor is the laser wavelength, 9 is the angle 
between the laser beams 19 , and <fi x .y,z is a phase detuning factor (both the latter 
are potentially variable). Such potentials have been realized in one 15 ' 16 , two (the 
so-called egg-carton potential) 20 and three dimensions 12 ' 21 . 

Moreover, present experimental realizations render feasible/controllable the adi- 
abatic or abrupt displacement of the magnetic or optical lattice trap 19 ' 22 ' 23 (induc- 
ing motion of the condensates), the "stirring" of the condensates providing angular 
momentum and creating excitations with topological charge such as vortices 24 ' 25 ' 26 
and vortex lattices thereof 27,28,29 . Additionally, phase engineering of the conden- 
sates is also feasible experimentally 30 , and this technique has been used in order to 
produce nonlinear matter-waves, such as dark solitons 31,32 in repulsive BECs. Note 
that, more recently, bright solitons have been generated as well 33,34 in attractive 
BECs, and both types are currently being studied extensively. 

Among these coherent structures, of particular interest are the nonlinear waves 
of non-vanishing vorticity. Vortices are worth studying not only due to their sig- 
nificance as a fundamental type of coherent nonlinear excitations but also because 
they play a dominant role in the breakdown of superfiow in Bose fluids 35,36,37 . The 
theoretical description of vortices in BECs can be carried out in a much more effi- 
cient way than in liquid He (see Ref. [38]) due to the weakness of the interactions 
in the former case. These advantages explain a large volume of work regarding the 
behavior of vortices in BECs, some of which has been summarized in Ref. [39]. 
It is interesting to note in this connection that the description of such topologi- 
cally charged nonlinear waves and their surprisingly ordered and robust lattices, 
as well as their role in phenomena as rich and profound as superconductivity and 
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superfluidity were connected to the theme of the recent Nobel prize in Physics in 
2003. 

It is around this exciting frontier of theoretical and experimental studies between 
atomic physics, soft condensed-matter physics and nonlinear dynamics that this 
brief review is going to revolve. Our aim is to report some recent developments 
in the study of vortices in the context of mean field theory (i.e., employing the 
GP equation). We consider various external potentials relevant to the trapping of 
the condensates, examining the existence, generation and dynamical stability of 
such coherent structures. It should be noted, however, that all the perturbations 
considered herein will preserve the Hamiltonian structure of the system. We will 
discuss these notions at various levels of increasing complexity, starting from that 
of a single vortex (in Section 3), proceeding to that of a few vortices and using 
the two- vortex system as a characteristic example (in Section 4), while in Section 
5, we will address the behavior of large clusters of vortices, the so-called vortex 
lattices. Section 6 summarizes our findings and presents some interesting directions 
for future studies. 



3. Single Vortex 
3.1. Generation 

It is well-known that if a superfluid is subjected to rotational motion, vortices will 
be generated in it. Such a situation also occurs in dilute BECs, where quantized 
vortices can be described in the framework of the GP equation. Thus, in this context, 
a vortex can typically be created upon "stirring" the condensate. In particular, 
beyond a critical angular velocity f2 c , the energy functional associated with the GP 
equation (incorporating the centrifugal term due to rotation), namely 



(4) 



is minimized by a single vortex configuration 40 , resulting in the generation of such 
a structure, observed also experimentally 25 . Note that in Eq. (4), 2 > f2 c is the 
angular velocity and L z — i(xd y — yd x ) represents the angular momentum along 
the z axis (* stands for complex conjugate). 

However, vortices can be spontaneously produced in a number of alternative 
ways, some of which we examine in what follows. More specifically: 

(1) they can spontaneously emerge as the pattern forming outcome of dynamical 
instabilities , or 

(2) they may be imprinted on the condensate via appropriate modulation of its 
phase. 

One instability that can be exploited as a method of producing vortex patterns 
is the so-called transverse or ("snaking") instability of rectilinear dark solitons. 
This instability, which has been studied extensively in the context of nonlinear 
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optics (see, e.g., Ref. [41, 42] for a review), forces a dark-soliton to undergo trans- 
verse modulations that cause the nodal plane to decay into vortex pairs. The snake 
instability is known to occur in trapped BECs as well (see Rcf. [32] for its experi- 
mental observation and Ref. [43] for relevant analytical and numerical results). In 
this context, dark solitons are placed on top of an inhomogeneous background, the 
so-called Thomas-Fermi (TF) cloud, which approximately yields the ground state 
wavefunction in the case of repulsive condensates (g > 0) 7 , and can be expressed 
as 



,t, i ■ ±\ /max{^- y ext ,0} 

*tf = exp(-z / ui)J , (5) 

where fi is the condensate's chemical potential. The snake instability of dark solitons 
in trapped BECs sets in whenever the soliton motion is subject to a strong coupling 
between the longitudinal and transverse degrees of freedom, i.e., far from ID ge- 
ometries. A relevant discussion demonstrating how the snaking instability manifests 
itself as the transverse confinement becomes weak, giving rise to the formation of 
vortices can be found in the recent work of [44]. 

To quantify better the above, we follow Ref. [45] to analyze the transverse 
instability via length-scale competition arguments. In particular, we first consider 
the following dimcnsionless 2D version (relevant for a quasi-two-dimensional, or 
"pancake" BEC lying on the x — y plane) of the original GP equation, 

** t = -^A ± * + |*| 2 * + K x tW*, (6) 

in which length is scaled in units of the fluid healing length £ = h/^/nogm (no 
is the peak density of the gas in the radial direction), t in units of £/c (where 
c = y/ nog/m is the Bogoliubov speed of sound) , and the atomic density is rescaled 
by the peak density n ; finally, the external potential is V ext {r) = (l/2)fi 2 r 2 , where 
r 2 = x 2 + y 2 and the parameter fl = y / u±Juj^(4Tral±no)~ 1 (where l± = y/h/mu)±_ 
is the transverse harmonic oscillator length, uj± being the transverse confining fre- 
quency) expresses the dimensionless effective magnetic trap strength. In the context 
of Eq. (6), and in the absence of the potential, the transverse instability occurs for 
perturbation wavenumbers 



k < C kr 



I \J sin 4 (f> + cos 2 <j) - (l + sin 2 (</>)) 



1/2 



(7) 



where cost/) is the dark-soliton amplitude (depth) and sm<fi is its velocity 41 . In the 
case of stationary (black) solitons cos cf> = 1, hence k cr — 1. On the other hand, in 
the presence of the potential, the characteristic length scale of the BEC (i.e., the 
diameter of the cloud in the TF approximation) is -Rbec = 2 3 / 2 ^o^ 2 /^, where fi is 
the dimcnsionless chemical potential. Then, one can argue that the criterion for the 
suppression of the transverse instability is that the scale of the BEC be shorter than 
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Fig. 1. Comparison of a subcritical case for xo = 10 with a supercritical one for xo = 15. The 2 
subplots on the left show different time snapshots of the evolution of the two-dimensional contour 
plot of the wavefunction square modulus for xq = 10. The 4 subplots of the right show the 
corresponding snapshots for xq = 15 in the unstable case, leading to the emission of a vortex pair. 



the minimal one necessary for the onset of the instability, leading to the condition 

n>^. (8) 

If inequality (8) is not satisfied, the transverse instability should develop, resulting 
in the breakup of the dipole configurations (resulting from a dark soliton, truncated 
by the Thomas-Fermi state) into vortex-antivortex pairs. It is relevant to note that 
similarly to the case of rectilinear solitons, the snake instability of the ring dark 
solitons can also be responsible for the creation of vortices and vortex arrays as well. 
In particular, as far as the ring dark solitons are concerned, they were previously 
predicted in the context of nonlinear optics 46 , where their properties were studied 
both theoretically 47 and experimentally 48 . These entities were also found to exist in 
the context of BECs 49 (see also the relevant recent work 50 ). In the latter context, 
and in the case of sufficiently large initial soliton amplitudes, ring dark solitons 
were observed to be dynamically unstable towards azimuthal perturbations that 
led to (snaking and) their breakup into vortex-antivortex pairs, as well as robust 
vortex arrays, the so-called "vortex necklaces" . The latter consist of four vortex- 
pair patterns, with their shape alternating between an "X" and a cross (for details, 
see Figs. 3 and 4 of Ref. [49]). Wc do not discuss these instabilities further, as they 
were analyzed in some detail in the recent review of Ref. [51]. 

A context similar to that of pattern forming instabilities, and one which may be 
loosely related to the Landau instability of superfiows, involves the interaction of the 
condensate with an impurity in a recently proposed dynamical experiment 52 (which 
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bears resemblances to the phenomenon of vortex trailing in fluids, see e.g., Ref. 
[53]). In particular, in Ref. [52] (and in the framework of the dimensionless GP Eq. 
(6)), the magnetic trap originally trapping the condensate at (0,0), was proposed 
to be displaced by a displacement of xq, but also in the presence of an additional 
anisotropic impurity potential of the form V^ mp = Vq sech 2 (^\/{x/r x ) 2 + (y/r y ) 2 ^ 
with Vq = 0.5, r y = 10 and r x — 5. If the displacement (which was giving rise to 
an oscillating motion of the condensate) was subcritical, then it was observed that 
it did not lead to the formation of "coherent structure radiation" . For supercrit- 
ical displacements, it was observed to give rise to the formation of vortex pairs, 
as shown in the results of Fig. 1. The critical speed (proportional to the critical 
displacement) acquired by the condensate upon "impact" with the impurity was 
numerically observed in Ref. [52] to closely match an effective speed of sound in the 
inhomogeneous medium (i.e., in the presence of the parabolic potential), indicating 
the possibility to attribute the vortex production in this context to a Landau-type 
instability 54 ' 55 . 




Fig. 2. Phase imprinting of vortices inside the combined magnetic and optical trap (left for vortices 
of charge 5 = 1) and solely in the magnetic trap (middle and right for 5 = 2 and 5 = 3 
respectively). The top and bottom row depict, respectively, the atomic density and phase. Left to 
right: vortices of increasing topological charge of one, two and three. In all cases, a radial magnetic 
trap frequency of fi = 0.1 has been used, while for the left panel the optical has been chosen with 
X x = X y = 4-7T, and t/> x = (fr y = 0, while Vb = 0.25. 

Let us conclude this section, by briefly discussing one of the standard techniques 
that have been implemented to create dark solitons and vortices in trapped BECs, 
namely the so-called phase imprinting or phase engineering technique. According to 
the latter, a dark soliton can be generated upon imprinting a phase difference of 7r 
along the condensate 30 ' 31 ' 32 . On the other hand, in the case of vortices, imprinting 
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(through an appropriate "phase mask" ) of a phase difference of 2tt around a contour 
can generate vortex structures (which carry topological charge). In fact, this method 
was proposed as a means of preparing vortex states in 2D BECs in Ref. [56] and 
was subsequently used in the laboratory in the experiments of Ref. [24] . Some of 
the advantages that this technique has over the previous ones are that: 

(1) It is very robust in the presence of even strong perturbations (i.e., it works 
equally well in the cases of combined potentials such as magnetic and optical 
ones), see left panels in Fig. 2. 

(2) It can be straightforwardly generalized to produce vortices of higher charge. 
These vortices may also be observed (for appropriate parameter values) to 
be dynamically stable for very long times and hence should, in principle, be 
experimentally observable; see middle and right panels in Fig. 2. 

3.2. Dynamical Stability 

3.2.1. Continuum Models 

Let us now address the question of stability of vortices trapped in a combined MT 
and OL as described in Ref. [57]. Consider a quasi- two-dimensional 12 condensate 
where a single vortex has been generated at the center of the combined potential 

V cxt (x,y) = Vmt(x,v) + V h(x,y), (9) 

where the contributions to the MT and OL, are given by the 2D equivalents of (2) 
and (3), respectively. As mentioned in the previous Section, the effective 2D GP 
equation [cf. Eq. (6)] applies to situations where the condensate has a nearly planar 
("pancake") shape, see for example Ref. [58] and references therein. Accordingly, 
vortex states considered below are not subject to 3D instabilities (corrugation of 
the vortex axis 39 ) as the transverse dimension is effectively suppressed. 

The stability of the vortex at the center of the trap can be qualitatively un- 
derstood in terms of an effective potential obtained by means of a variational 
approximation 59 . As in Ref. [60], we use the following ansatz to approximate the 
position ro(t) = (xo(t), yo(t)) of vortex near the trap center 

*(x,y,t) = B(t)\\ro(t)\\exp[-\\r (t)\\ 2 /b(t)] e^°« (10) 

where (fio(t) = tan _1 [(y — y (t))/(x — xo(t))] and ||ro(t)|| denotes the 2-norm of the 
vector tq. Similarly to the calculations in Rcfs. [61, 62], or upon employing the con- 
servation of norm, it is straightforward to show that, to leading order, B{t) = B(0) 
and b(t) = 6(0) are approximately constant. Then, assuming the same detuning 
and periodicity in all directions (i.e., <f> x = (ft y = (ft and X x = X y = A), the substitu- 
tion of the ansatz (10) into the Lagrangian form of the GP equation (6) leads to a 
quasi-particlc description of the vortex center through the Newton-type equations 
of motion 

dV e s(x ,yo) , .. dV e s(x ,y ) 

x = j , and y Q = , (11 

dx dy 



I 



I 
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X X »io" s X X 



Fig. 3. Vortex stability inside a combined magnetic and optical trap. The MT and OL parameters 
are u x = ujy = 0.002, Vb = 0.5 and A = \ x = \y = 2n. Left: the vortex is stable at the bottom 
of a cosinusoidal OL [<f> = in Eq. (3)]. The top panel shows the contour plot of the density at 
t = 100 (138 ms). The bottom left panel is a cut of the same density profile along x = 0, while 
the bottom right one shows the motion of the vortex center for < t < 100, the initial position 
being marked by a star. Notice the scale (10~ 3 , or 1 nm in physical units) of very weak motion 
of the vortex, which thus stays practically immobile at the origin. Right: same as left panels but 
with a sinusoidal OL [<j> = tt/2] and for a larger time of t = 250 (346 ms). The bottom right panel 
depicts the motion of the vortex center for < t < 250 [positions of the vortex center at t = 100 
(138 ms) and t = 200 (276 ms), respectively, are indicated by the star and circle]. 

where the effective potential is given by 

+ i^(x 2 + 2/ 2 ), (12) 

and Q((f>) is given by a rather cumbersome expression. In our case of interest we 
will only need the following values 

0(0) = £ - l) exp {^f) , and Q (1) = -Q(0), (13) 

Equations (11) and (12) indicate that the coordinates of the vortex center are 
prescribed by two uncoupled nonlinear oscillators (see also Ref. [63] for a similar 
result in the context of optics). The above, generalizes the well-known result 64 
that the center of a dark soliton (the ID "sibling" of the vortex) behaves like a 
Newtonian particle in the presence of the external potential (see relevant work for 
dark matter- wave solitons in optical lattices in Ref. [65]). Figure 3 depicts the vortex 
stability for the two detuning cases in (13). For A = 2n, Vb = 0.5 > and 6=1 
(by fitting ansatz to numerical solution), we have that Q(<f> = 0) = — Vq/8 < and 
Q{<j> = tt/2) = —Q(0) = Vq/8 > 0. Therefore, a cosinusoidal OL ((/> — 0) produces 
a stable vortex (see left panels of Fig. 3) while a sinusoidal OL (tp — tt/2) induces 
an instability (see right panels of Fig. 3). Note that, recently, relevant results have 
been obtained using different approaches 66 . 

It is worth mentioning that the variational approach outlined above, although 
capable of capturing the stability of the vortex at the center of the trap, fails to 



V eB {x,y)=Q(<l>) 



( Attx \ ( Airy 

cosf — )+cos( — 
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reproduce the correct dynamics. More specifically, as shown in Fig. 3, the vortex 
center follows an outward spiraling motion. This spiral motion is the combination of 
the unstable behavior captured by the Hamiltonian dynamics of Eqs. (11) and (12) 
together with the well-known precession of vortices inside the MT 67 . The precession, 
not captured by our approximation, predicts that, close to the trap center, a vortex 
rotates with an angular frequency given by 67 

-3u x u y ( Rbec\ n .s 
^prcc = — ^— - In I — — I , (14) 

where (i is the chemical potential, i?BEC is the mean transverse dimension of the 
condensate and £ is the width of the vortex core, which, in fact, is the same as the 
healing length of the condensate 7 . 

3.2.2. Discrete Models 

An interesting situation arises for strong optical lattices (Vo ^ /j) where the con- 
densate is effectively transformed into a collection of "droplets" (in each lattice well) 
that can be described by the spatially discrete analogue of the NLS 15 ' 17,18 ' 21,68 . In 
this case, it is possible to describe the evolution of the condensate wave function 
by the discrete nonlinear Schrodinger equation (DNLS) as can be seen through a 
Wannier function decomposition 69 ' 70 . In non-dimensional units, the DNLS reads 

l^ + CA^cfln + ttnUn^O, (15) 

where <j). q is the condensate wave function localized at the OL trough with coor- 
dinates rj — (m,n) and 77 = (l,m,n), for the 2D and 3D cases respectively, C 

is the coupling constant, and A^ stands for the discrete Laplacian in d diluen- 
ts (3) 

Sions: '(j) m ,n = <f>m+l,n + <t>m.n+l + <frm,n-l + <Pm-l,n - Hm.n and A 2 </>;, TO ,„ = 
4>l+l,m,n+<j>l,m+l,n+<j>l,m,n+l+4 l l~l,m,n + 4>l,m~l,n-l+4>Lm,n-l— §4>l,m,n- In Rcfs. [71, 

72, 73] stationary solutions of (15) are sought by considering (f> v = cxp(— ifiot) u v , 
where /xo is the dimensionless chemical potential of the condensate, that yields to 
the time-independent equation: 

-^0«r, = CA^Ur, + \u, q \ 2 U n , (16) 

where \u v \ 2 is proportional to the atomic density at the 77-th trough. Since Eq, 
(16) has a scale invariance, ft can be fixed arbitrarily. As in Refs. [71, 72, 73], by 
using an appropriate initial discrete vortex ansatz, it is possible to find numerical 
solutions to Eq. (16). Then, by applying continuation-type methods, together with 
linear stability computations, the branches of existence and stability for discrete 
vortices of different charges in two- (d = 2) and three-dimensional (d = 3) settings 
can be obtained. 

The stability results presented in Rcfs. [71, 73], for vortices of charge S (e.g., 
for dimensionless chemical potential /io = —4), can be summarized in the following 
stability table: 
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d 


5 = 


5 = 1 


5 = 3 


5 = 2 


2 


C < C c ( r 2 ' 0) « 4 


C < C^ ,1} « 1.6 


C < d?' 3) » 0.398 


U 



where the approximate region for stability is indicated (the letter "U" denotes 
unstable for all values of C). As indicated, no stable vortices with 5 = 2 were 
obtained 74 . Nonetheless, Eq. (16) admits real stationary solutions — generated by 
the real part of the corresponding genuine vortex (which is complex) — that are sta- 
ble for sufficiently weak coupling. These real solutions, the so-called quasi- vortices, 
correspond to quadrupoles (5 = 2) and octupoles (5 = 4) 71,73 . Similar results were 
found also for the three-dimensional case in Ref. [72]. In Fig. 4 we depict a few 
examples of discrete vortices for different vorticities in two- and three-dimensions. 




Fig. 4. Vortices in the discrete nonlinear Schrodinger equation in two (d = 2) and three (d = 3) 
dimensions. The top (bottom) row depicts the real (imaginary) part of the stationary solution (cf. 
Eq. (16)) with vorticity S. For d = 3 (b,c,d) the plots depict the contour plots for Re(if( m ,n) = 
iconstant (top row) and Im(«; m n ) = iconstant (bottom row). The parameters are, from left 
to right: (a) d = 2, S = 3, C = 0.02 (b) d = 3, S = 1, C = 0.7, constant=0.5 (c) d = 3, S = 2, 
C = 0.01, constant=0.25 and (d) d = 3, S = 3, C = 0.01, constant=0.25. 



It is important to mention that the discreteness is responsible for inducing the 
stability of vortices in three-dimensional settings that otherwise (in the continuum 
model) are strongly unstable. A natural question that arises is the fate of unstable 
solutions. For example, we have observed that, in two dimensions, a vortex with 
5 = 3 and C = 0.618 > Cg' may decay into a configuration consisting of a 
combination of a stable soliton (5 = 0) and a stable 5 = 1 vortex (here C < ci? 
and C < Ccr' )■ On the other hand, we have observed a striking effect where, 
in three-dimensions, an unstable 5 = 2 vortex decays into a stable 5 = 3 vortex 
(for C = 0.01 < Cl 3 ' 3) ), thereby increasing the total topological charge instead of 
decaying to a combination of lower order (5 = 0, 1) vortices. It should be noted that 
this change of vorticity is possible in the discrete lattice model, in which the angular 
momentum is not a dynamical invariant. Another noteworthy vortex solution also 
found in Refs. [72] consists of a complex of two mutually orthogonal 5=1 vortices 
(one in each component) in the discrete version of a two-component condensate (cf. 
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Section 4.2). 

4. Two Vortices 

4.1. One- Component BEC: Interactions 

Up to this point we have dealt with the generation and stability of vortices. Let 
us now focus on the important issue of vortex-vortex interactions. As is the case 
also for fluid vortices, two BEC vortices of opposite charge travel parallel to each 
other at constant speed c, while vortices of same charge rotate at constant angular 
speed a 75 ' 76 . From direct numerical integration of Eq. (6) we have been able to 
characterize the dynamics of interacting vortices in the absence of any trapping 
potential (fi = 0). With an appropriate initial phase mask (see Section 3.1) we seed 
two vortices, with respective vorticities S\ and S2, at a desired distance from each 
other. Then, we fitted a Pade approximation 77 to find the vortex centers during 
evolution, and numerically extracted the linear (angular) velocities for vortices of 
opposite (equal) charge as a function of their separation distance p (see left panel in 
Fig. 5). We performed our experiments for same and opposite charge vortices with 
vorticities |Si| = IS2I = 1 and \S±\ = IS2I = 2. The angular and linear velocities of 
interacting vortices seem to nicely obey power laws for a wide range of separations 
(10 < p < 50): c(p) ~ p~ x and a(p) ~ p~ 2 (see left panel of Fig. 5). 




In(p) 



Fig. 5. Left: linear c(p) (angular ct(p)) velocity for two vortices of opposite (equal) charge |5| = 
I Si I = IS I = 1 and \S\ = I Si I = IS2I = 2 as a function of their separation distance p. Middle and 
right panels depict two configurations of three interacting vortices of charge S = 1. The solid and 
dashed line represent the trajectories of their centers from direct numerical integration of (1) and 
from the quasi-particle approximation of (18), respectively. 

The next step to describe the dynamics of interacting vortices consists on iden- 
tifying an appropriate Lagrangian that gives rise to the correct dynamics. In order 
to achieve this, it is crucial to note that BEC vortices are known to obey first or- 
der, kinematic equations 39 , contrary to what is known for the Newtonian (second 
order) equations describing other types of solitary wave center of mass dynamics 59 . 
To circumvent this obstacle, wc need to construct an interaction Lagrangian that 
gives the "correct" vortex dynamics through its Euler-Lagrange equations. Such a 
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Lagrangian for two vortices, with same vorticity S = S m = S n , may be given in 
the form 

L m ,n ~ det(r„, r m ) +det(r TO , r„) - ASln(p) (17) 

where f m = {x m ,y m ) T is the position vector of vortex m and A is a constant. For 
vortices of opposite charge one needs to multiply the velocities for each vortex by 
their respective charge. This is equivalent to including the vortex charge in the 
definition of the Poisson brackets for the vortex interaction Hamiltonian (cf., Ref. 
78 and references therein). 

Let us explicitly write the equations of motion for same charge vortices. In this 
case, the pairwise Lagrangian (17) gives rise to the well known equations of motion 
for two (fluid) point vortices centered at f n = (x n ,y n ) T and r m — (x m ,y m ) T : 

■ __«n Dm ~ Dn 
X m - AQ 2f) 2 

i A O X ' m ~ Xn 

y m = +Ab 



2p 



2 



where p = (x rn — x n ) 2 + (y m — y n ) 2 is inter-vortex distance and A — A\ 3.96 
for S = 1 and A = A 2 « 7.80 for S — 2 is a constant determined from the numerics 
(cf. left panel in Fig. 5). As expected we have 2Ai = 7.92 s=a A 2 , namely, the effects 
of an 5* = 2 vortex are equivalent to the superposition of the effects of two nearby 
S = 1 vortices. For opposite charge vortices traveling parallel to each other at 
constant speed c, our numerics predict that c(p) = Bp -1 with B = B\ s=s 2.15 
and B = B 2 ~ 4.43 for singly- (l^l = 1) and doubly-charged (l^l = 2) vortices 
respectively. As for same charge vortices, the relationship 2B\ = 4.30 « B 2 holds 
approximately. 

Approximation (18) treats vortex pairs as quasi-particles with interacting po- 
tentials. It is important to mention that this approximation to the dynamics dete- 
riorates when the vortices get too close to each other. Nonetheless, we have checked 
that the approximation remains valid down to separation distances of the order of 
the healing length (i.e., the width of the vortices) for | | = 1 vortices. For |,S| = 2, 
the two vortices tend to split into two pairs of £ = 1 vortices 79 when the separation 
distance was decreased below approximatively p = 10. 

In the middle and right panel of Fig. 5 we present two examples with three 
interacting vortices. The solid line represents the actual orbits obtained from direct 
numerical integration of Eq. (1), while the dashed line depicts the orbits using 
the quasi-particle approximation (18). As it can be observed from the figure, the 
superposition of all the interactions, given by (18), for the three vortices is in good 
agreement with the full model. It is worth mentioning that the case depicted in the 
right panel of the figure corresponds to an orbit with an initial condition close to 
an unstable configuration. This explains the larger discrepancy for this case when 
compared to the middle panel. 

We note in passing that a more general approach towards computing vortex 
interactions may involve the use of a Lagrange multiplier at the level of the PDE 
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(initialized with two vortices) which maintains the distance between the vortices 
fixed. Then, the force associated with the relevant multiplier is what is needed to 
balance the interaction force between the vortices and hence can also be used to 
infer the interaction potential. 

4.2. Two-Component BECs: Dipole Bound States 

A very relevant generalization of the class of physical systems that we have dis- 
cussed so far, and of the solitary waves they can support, concerns the case of 
coupled multi-component BECs. There has recently been a considerable volume of 
work relevant to the properties of coupled BECs ranging from the study of ground 
state solutions 80 ' 81 to small-amplitude excitations 82 . Furthermore, the formation of 
various structures such as domain-walls 83 ' 84 ' 85,86 , bound dark-dark 86 , dark-bright 87 , 
dark-antidark, dark-gray, bright-antidark and bright-gray soliton complexes 88 , as 
well as spatially periodic states 89 was also predicted. On the other hand, experi- 
mental results have been reported for mixtures of different spin states of 87 Rb (see 
Ref. [90]) and mixed condensates 91 ' 92 . It is relevant to also mention the efforts to- 
wards the realization of two-component BECs from different atomic species, such 
as 41 K- 87 Rb (see Ref. [93]) and 7 Li- 133 Cs (see Ref. [94]). 

Typically, the generalized mean field model for two coupled BECs involves two 
nonlinearly coupled CP equations. However, in experiments with a radio-frequency 
(or an electric field) coupling two separate hyperfine states 90 , the relevant model 
also involves a linear coupling between the wavefunctions. The governing normalized 
equations are then of the form: 

1 



iip2t 



-A + V + g 11 \tp 1 \ 2 +g 12 \^ 2 \ 2 
-^A + V + g 12 \tpi\ 2 + g 2 2\ip2\ 2 



fa+anfa, (19) 



ih + on/>i, (20) 



where the intra- and inter-species interactions are characterized by the coefficients 
9jj (j ~ 1)2) and gi 2 respectively, while a denotes the strength of the radio- 
frequency (or electric field) coupling. In the recent work of Ref. [95], the special 
case of <7n = 522 = .912 = g was examined. The latter can be written in a vector 
form as: 

iipt - aP4> = -^AiP + {^G^)tp + V{x)ij, (21) 



where 



1 

1 



(22) 



and G = gl, with I being the unit matrix. In that case, as was also previously 
known in optics 96 , one can make a unitary transformation: 

/ tt^\ 1 -iaPt / ( cos(crf) —ism(at)\ , . . 

\— ism(atj cos(at) / 



I 



I 
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Fig. 6. Left panels: Contour plots of \tpi\ 2 for two coupled vortices, initially placed at x = ±5, 
for t = (a), T/4 (b), T/2 (c), and 3T/4 (d), with T = 7r/a Ri 15.7 (a = 0.2); H = 0.045, 
A = 5HS22 — = — 9 x 10 -4 ( 87 Rb). The vortices "interchange locations" (in a structure 
resembling a spiral wave). Right panels: Same as the left but for A = — 3 (g = 1, h = 2). The 
configuration breaks up forming spiral patterns. 



Then the original set of equations is transformed into: 

i^ot - -\Mo + (V'jGWo + ^(x)Vo, (24) 

i.e., the linear coupling can be completely eliminated from the equations. Notice 
that this special case of approximately equal scattering length coefficients is relevant 
to the experiments performed with two spin states of 87 Rb (see Ref. [91]), where 
5ii : ffi2 : 522 = 1-03 : 1 : 0.97. 

In the system of nonlinearly coupled GP equations for ?/> , one can construct a 
dipolc configuration with a pair of vortex structures, see Fig. 6. This configuration 
was obtained in the figure, by means of imaginary time integration in the absence 
of linear coupling, starting with one component having a vortex centered at (5, 0), 
while the other has a vortex at (—5, 0). After the configuration relaxes to the sta- 
tionary vortex pair solution of Eqs. (19)-(20), one can turn on the linear coupling 
and obtain a spiral rotation between the vortices resembling a spiral wave. While 
this solution is exact for h = 312/511 = 1 (assuming 311 = 322), it also persists for 
non-unit values of h. In particular, the spiral structure persists for h < h c = 1.32 
beyond which the regularity of the Rabi oscillations of matter between the compo- 
nents is destroyed. In this case, the breakup leads to the formation of spiral patterns 
in the condensate 95 . 

5. Vortex Lattices 

We now turn to a "thermodynamic limit" level which is also, however, particularly 
relevant to BEC experiments. Rapid rotation of 2D BECs induces the generation 
of many vortices 39 ' 97 ' 98 that typically settle into ordered lattices 25 - 26 ' 27 ' 28 ' 29 ' 99 ' 100 . 
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Particularly enticing in this respect are the available pictures of such lattices and 
their (practically perfect) triangular patterns 101 . These are the so-called Abrikosov 
lattices 102 , that were long ago predicted in the theory of superconductivity (and 
that are cited as the prediction that earned their discoverer the Nobel prize in 
Physics in 2003). In the context of type-II superconductors, free energy arguments 
can be used to demonstrate that the triangular lattice is the most energetically 
favorable (ground-state) configuration 103 . 

It is quite interesting in this context, to study vortex lattices in the framework 
of Eq. (6) both in the presence of the external magnetic trap (MT), as well as in the 
one of the optical lattice (OL) . Naturally, it is relevant to perform direct simulations 
of such lattices for different (external potential) parameter values. However, a quite 
relevant alternative tool in order to study the ground states of such configurations 
and their structural transitions is the use of molecular dynamics (MD) techniques 
such as the Parrinello- Rahman (PR) method 104 . In the PR dynamics, not only are 
the positions of the vortices (viz. particles) evolved in time, but so are the coordinate 
vectors of the box in which the coherent structures are located. In particular, if the 
system of coordinates consists of the vectors a = (a x , a v ) and b = (b x , b y ), then the 
metric tensor becomes (in 2D) G = h T h, where h is the coordinate transformation 
matrix with a, b as its rows [(x m , y m ) T = h-(£ n , ?7n) T ]- Then the Lagrangian for the 
fictitious coupled dynamics of the "particle" -lattice and the MD box reads 104,105 : 



where M,W are the particle and box mass respectively, (£„,f„) are the coordi- 
nates of the n-th particle in the "box frame", V mn is the interaction potential 
between quasi-particles m and n, and the summations are over the total number 
of vortices N v . Then, one can perform PR-MD by solving the ensuing dynami- 
cal equations 104,105 , to obtain the coordinate system evolution (along with that of 
the "particle" positions). More specifically, from (25) one obtains the dynamics for 
each of the AN V + 8 degrees of freedom, say q n and q n , via the Euler-Lagrangc 
equations d/dt(dC/dq n ) = dC/dq n . In order to emulate the thermodynamic limit 
of the system (large number of vortices), the minimal image convention can be 
employed where the box is periodically repeated for all of its 8 neighbors and, for 
each vortex, the N v largest contributions, from all the 9 boxes, are taken into ac- 
count. The results obtained in Ref. [106] indicate that, interacting spiral vortices in 
a Ginzburg-Landau field-theoretic context 107 (which, in fact, is a dissipative version 
of the GP equation — see also Ref. [108]), typically settle, for sufficiently large num- 
ber of vortices, to configurations similar to the experimentally observed triangular 
configuration (sec Fig. 7). Note that similar results can be applied to pulses where 
the topological charge is zero. 

In order to specifically apply the PR-MD simulation to a gas of vortices in the 




n 



(25) 



+ -W(dl+d 2 y +bl+b 2 y )-Y,Vr 



i 
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Fig. 7. Typical vortex lattice MD simulation with a random initial pattern of N v = 225 same- 
charge vortices. After a short transient (t < 8), the configuration equilibrates into a triangular 
lattice (t = 10). 



BEC model (1), one needs to recast the MD Lagrangian (25) using the realistic 
vortex-vortex Lagrangian (17) 109 . The form of the Lagrangian (25) allows for a 
direct incorporation of the effects induced by the MT and OL. It is also possible to 
study the configuration changes in the presence of quadrupolar excitations or other 
symmetry stresses exerted on the vortex lattice (motivated by the experiments 
of Refs. [29, 99, 100]). We should note that the external imposition, in the box 
containing the vortices, of shear stresses of different symmetry is quite feasible 
experimentally. It can be implemented through a rotating laser manipulating the 
boundary of the box (in this case, the periphery of the Thomas-Fermi cloud) 110 . 




Fig. 8. Molecular dynamics simulations with an initial pattern (a) corresponding to a slightly 
perturbed square lattice of N v = 225 same-charge vortices, (b) After some transient, the config- 
uration equilibrates to a star pattern, (c) Rhomboidal pattern obtained from a slightly different 
initial perturbation of (a). 



The PR-MD approach is therefore a powerful tool that can be used to identify 
the conditions for which the triangular state persists as the ground state configura- 
tion (for such conditions, see for example Ref. [111]). It can also be used to obtain 
the structural transition points (obtained through direct simulations, e.g., in the 
presence of the optical lattice in Ref. [112]) to other states such as rhomboidal or 
star-like patterns (see for example Fig. 8) . The approach consists of analyzing the 
dynamically obtained PR-MD states through stability computations, standard con- 
tinuation and bifurcation theory tools to explore the structural phase transitions 
and to obtain the instability eigenvectors that lead to phase transformations. An al- 
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ternative dynamical scenario, very relevant to recent experimental settings, involves 
the annihilation of a central chunk of vortex lattice matter, through a localized 
laser heating 100 , that results in the remaining vortex lattice exhibiting oscillatory 
modes, known as Tkachenko oscillations. Such modes may also be identified 109 as 
limit-cycle, time-periodic solutions of the PR-MD numerical procedure. 

6. Summary and outlook 

While there is a large volume of work on the waves of topological charge in Bosc- 
Einstcin condensates (a large fraction of which is summarized in Ref. [39] , as well 
as in the present brief review), there are still numerous open problems regarding 
the vortex state that this experimentally and theoretically tractable context may 
allow us to explore. 

Clearly, a prominent position among such open problems (in the context of 
isolated vortices) is the question of a mathematical understanding of the detailed 
dynamical stability picture of vortices of various topological charges (S > 1) in 
the presence of the combined magnetic and optical trappings. A first step in that 
direction is offered by the work of Rcf. [113]; however there are still many open 
questions concerning the effect of the trapping potentials. 

Another subject that apparently has received very little attention and whose 
theoretical understanding is still to a large extent incomplete (in the context of few 
vortices) is the behavior of vortex dipoles. Such configurations have been now ob- 
tained for two-component BECs in the work of Refs. [95, 114], however topics such 
as the interaction of such dipoles and the ensuing dynamics are still unexplored. 
From the mathematical point of view, a first attempt at obtaining reduced equa- 
tions that adequately describe dipole dynamics has been given in Ref. [115]. The 
applicability of such an approach in the context of BECs would be of particular 
interest. 

There are also numerous interesting open questions regarding vortex lattices. 
The dynamical imposition/time dependence of external potentials in conjunction 
with possible interaction between multiple condensate components may produce 
interesting effects of frustration that may induce structural phase transitions such 
as the ones discussed in Ref. [112] (e.g., a reshaping from a triangular lattice to 
a square one; see for example Ref. [116]) The PR-MD setup is straightforwardly 
amenable to the inclusion of such external effects. By properly incorporating a 2D 
OL, we expect to induce a locally attractive (or locally repulsive) energy landscape 
for each vortex at any required location. This should enable us to engineer a rich va- 
riety of "target" lattice configurations — provided that their energy is not far from 
a local minimum. Applications of this (as well as statistical mechanics) techniques 
to understand the ground state of vortex lattices under external perturbations or 
multi-component interactions is bound to provide interesting conclusions for the 
non-equilibrium thermodynamics of such multi-vortex patterns. Notice that some 
of these tasks (such as identifying stationary vortex lattice states and computing 
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the corresponding eigenvalues of linearization around them) can be performed by 
methods that have been developed in matrix-free numerical linear algebra 117 . This 
can be done by using appropriately initialized short bursts of time evolution simu- 
lations instead of the very expensive large Jacobian eigenvalue computations. 

These are only some among the many questions/topics that are now starting to 
be addressed (for instance, one can ask the same questions at finite temperature 
and try to understand the interaction of the vortex condensate with the gas in that 
case). We hope that this intriguing journey still hides, as Cavafy says in his Ithaca, 
a lot of "ports seen for the first time" . . . 
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